home *** CD-ROM | disk | FTP | other *** search
/ CU Amiga Super CD-ROM 19 / CU Amiga Magazine's Super CD-ROM 19 (1998)(EMAP Images)(GB)[!][issue 1998-02].iso / CUCD / Programming / LEDA / prog / plane / seg_int.c < prev    next >
C/C++ Source or Header  |  1994-08-05  |  15KB  |  633 lines
  1. #include <LEDA/sweep_segments.h>
  2.  
  3. #include <LEDA/segment.h>
  4. #include <LEDA/graph.h>
  5. #include <LEDA/prio.h>
  6.  
  7. #include <LEDA/sortseq.h>
  8. #include <LEDA/impl/skiplist.h>
  9. #include <LEDA/impl/ab_tree.h>
  10.  
  11. #include <math.h>
  12.  
  13.  
  14. #define EPS  0.00001
  15. #define EPS2 0.0000000001
  16.  
  17. class sw_POINT;
  18. class sw_SEGMENT;
  19. typedef sw_POINT* sw_point;
  20. typedef sw_SEGMENT* sw_segment;
  21.  
  22. enum sw_point_type {Cross=0,Rightend=1,Leftend=2}; 
  23.  
  24. class sw_POINT
  25. {
  26.   friend class sw_SEGMENT;
  27.  
  28.   sw_segment seg;
  29.   int     kind;
  30.   double  x;
  31.   double  y;
  32.  
  33.   public:
  34.  
  35.   sw_POINT(double a,double b)  
  36.   { 
  37.     x=a; y=b; seg=0; kind=Cross;
  38.    }
  39.  
  40.  
  41.  
  42.   sw_POINT(point p)         
  43.   { 
  44.     x=p.xcoord();y=p.ycoord();seg=0;kind=Cross;
  45.    }
  46.  
  47.  
  48.   LEDA_MEMORY(sw_POINT);
  49.  
  50.   friend double     get_x(sw_point p);
  51.   friend double     get_y(sw_point p);
  52.   friend int        get_kind(sw_point p);
  53.   friend sw_segment get_seg(sw_point p);
  54.  
  55.   friend bool intersection(sw_segment, sw_segment, sw_point&);
  56. };
  57.  
  58. inline double     get_x(sw_point p)    { return p->x; }
  59. inline double     get_y(sw_point p)    { return p->y; }
  60. inline int        get_kind(sw_point p) { return p->kind; }
  61. inline sw_segment get_seg(sw_point p)  { return p->seg; }   
  62.  
  63.  
  64.  
  65. static int compare(sw_point& p1,sw_point& p2)
  66. { if (p1==p2) return 0;
  67.  
  68.   double diffx = get_x(p1) - get_x(p2);
  69.   if (fabs(diffx) > EPS2 ) return (diffx > 0.0) ? 1 : -1;
  70.  
  71.   int    diffk = get_kind(p1)-get_kind(p2);
  72.   if (diffk != 0) return diffk;
  73.  
  74.   double diffy = get_y(p1) - get_y(p2);
  75.   if (fabs(diffy) > EPS2 ) return (diffy > 0.0) ? 1 : -1;
  76.  
  77.   return 0;
  78.  
  79. }
  80.  
  81.  
  82. class sw_SEGMENT
  83. {
  84.   sw_point startpoint;
  85.   sw_point endpoint;
  86.   double  slope;
  87.   double  yshift;
  88.   node  left_node;
  89.   int   orient;
  90.   int   color;
  91.   int   name;
  92.  
  93.  
  94.   public:
  95.  
  96.   sw_SEGMENT(sw_point, sw_point,int,int);     
  97.  
  98.  ~sw_SEGMENT() { delete startpoint; delete endpoint; }     
  99.  
  100.   LEDA_MEMORY(sw_SEGMENT);
  101.  
  102.   friend sw_point get_startpoint(sw_segment seg);
  103.   friend sw_point get_endpoint(sw_segment seg);
  104.   friend double get_slope(sw_segment seg);
  105.   friend double get_yshift(sw_segment seg);
  106.   friend int  get_orient(sw_segment seg);
  107.   friend int  get_color(sw_segment seg);
  108.   friend int  get_name(sw_segment seg);
  109.   friend node get_left_node(sw_segment seg);
  110.   friend void set_left_node(sw_segment seg, node v);
  111.  
  112.   friend bool intersection(sw_segment, sw_segment, sw_point&);
  113. };
  114.  
  115. inline sw_point get_startpoint(sw_segment seg)    { return seg->startpoint; }
  116. inline sw_point get_endpoint(sw_segment seg)      { return seg->endpoint; }
  117. inline double get_slope(sw_segment seg)           { return seg->slope; }
  118. inline double get_yshift(sw_segment seg)          { return seg->yshift; }
  119. inline int  get_orient(sw_segment seg)            { return seg->orient; }
  120. inline int  get_color(sw_segment seg)             { return seg->color; }
  121. inline int  get_name(sw_segment seg)              { return seg->name; }
  122. inline node get_left_node(sw_segment seg)         { return seg->left_node; }
  123. inline void set_left_node(sw_segment seg, node v) { seg->left_node = v; }
  124.  
  125.  
  126. sw_SEGMENT::sw_SEGMENT(sw_point p1,sw_point p2,int c, int n)    
  127.   {
  128.     left_node  = nil;
  129.     color      = c;
  130.     name       = n;
  131.  
  132.     if (compare(p1,p2) < 0)
  133.      { startpoint = p1; 
  134.        endpoint = p2; 
  135.        orient = 0;
  136.       }
  137.     else
  138.      { startpoint = p2; 
  139.        endpoint = p1; 
  140.        orient = 1;
  141.       }
  142.  
  143.     startpoint->kind = Leftend; 
  144.     endpoint->kind = Rightend; 
  145.     startpoint->seg = this; 
  146.     endpoint->seg = this;
  147.  
  148.     if (endpoint->x != startpoint->x)
  149.     {
  150.       slope = (endpoint->y - startpoint->y)/(endpoint->x - startpoint->x);
  151.       yshift = startpoint->y - slope * startpoint->x;
  152.  
  153.       startpoint->x -= EPS;
  154.       startpoint->y -= EPS * slope;
  155.       endpoint->x += EPS;
  156.       endpoint->y += EPS * slope;
  157.     }
  158.     else //vertical segment
  159.     { startpoint->y -= EPS;
  160.       endpoint->y   += EPS;
  161.      }
  162.   }
  163.  
  164.  
  165. static double x_sweep;
  166. static double y_sweep;
  167.  
  168.  
  169. static int compare(sw_segment& s1,sw_segment& s2)
  170. {
  171.   double y1 = get_slope(s1)*x_sweep+get_yshift(s1);
  172.   double y2 = get_slope(s2)*x_sweep+get_yshift(s2);
  173.  
  174.   double diff = y1-y2;
  175.  
  176.   if (fabs(diff) > EPS2) return (diff > 0.0) ? 1 : -1;
  177.  
  178.   if (get_slope(s1) == get_slope(s2)) 
  179.         return compare(get_x(get_startpoint(s1)), get_x(get_startpoint(s2)));
  180.  
  181.   if (y1 <= y_sweep+EPS2)
  182.         return compare(get_slope(s1),get_slope(s2));
  183.   else
  184.         return compare(get_slope(s2),get_slope(s1));
  185.  
  186. }
  187.  
  188.  
  189. void Print(sw_segment& x) 
  190. { sw_point s = get_startpoint(x);
  191.   sw_point e = get_endpoint(x);
  192.   cout << 
  193.     string("[(%f,%f)----(%f,%f)]",get_x(s),get_y(s), get_x(e),get_y(e));
  194. }
  195.  
  196.  
  197. static priority_queue<seq_item,sw_point>* X_structure;
  198.  
  199. static sortseq<sw_segment,pq_item>* Y_structure;
  200.  
  201.  
  202.  
  203. bool intersection(sw_segment seg1,sw_segment seg2, sw_point& inter)
  204. {
  205.   if (seg1->slope == seg2->slope)
  206.     return false;
  207.   else
  208.   {
  209.     //x-coordinate  of the intersection
  210.     double Cross_x = (seg2->yshift - seg1->yshift) / (seg1->slope - seg2->slope);
  211.  
  212.     if (Cross_x <= x_sweep) return false;
  213.  
  214.     double s1 = seg1->startpoint->x;
  215.     double s2 = seg2->startpoint->x;
  216.     double e1 = seg1->endpoint->x;
  217.     double e2 = seg2->endpoint->x;
  218.  
  219.     if (s2>Cross_x || s1>Cross_x || Cross_x>e1 || Cross_x>e2) return false;
  220.  
  221.     //y-coordinate of the intersection
  222.     double Cross_y = (seg1->slope * Cross_x + seg1->yshift);
  223.     inter =  new sw_POINT(Cross_x,Cross_y);
  224.  
  225.     return true;
  226.   }
  227. }
  228.  
  229.  
  230. static pq_item Xinsert(seq_item i, sw_point p) 
  231. { 
  232.   return X_structure->insert(i,p);
  233. }
  234.  
  235. static sw_point Xdelete(pq_item i) 
  236. {
  237.   sw_point p = X_structure->inf(i);
  238.   X_structure->del_item(i);
  239.   return p;
  240. }
  241.  
  242.  
  243. static node New_Node(GRAPH<point,int>& G,double x, double y )
  244. { return G.new_node(point(x,y)); }
  245.  
  246. static void New_Edge(GRAPH<point,int>& G,node v, node w, sw_segment l )
  247. { if (get_orient(l)==0)
  248.        G.new_edge(v,w,get_color(l));
  249.   else G.new_edge(w,v,get_color(l));
  250. }
  251.  
  252.  
  253. void vertical_segment(GRAPH<point,int>& SUB, sw_segment l)
  254. { 
  255.   sw_point p = new sw_POINT(get_x(get_startpoint(l)),get_y(get_startpoint(l)));
  256.   sw_point q = new sw_POINT(get_x(get_endpoint(l)),get_y(get_endpoint(l)));
  257.  
  258.   sw_point r = new sw_POINT(get_x(p)+1,get_y(p));
  259.   sw_point s = new sw_POINT(get_x(q)+1,get_y(q));
  260.  
  261.   sw_segment bot = new sw_SEGMENT(p,r,0,0);
  262.   sw_segment top = new sw_SEGMENT(q,s,0,0);
  263.  
  264.   seq_item bot_it = Y_structure->insert(bot,nil);
  265.   seq_item top_it = Y_structure->insert(top,nil);
  266.   seq_item sit;
  267.  
  268.   node u,v,w;
  269.   sw_segment seg;
  270.   
  271.  
  272.   for(sit=Y_structure->succ(bot_it); sit != top_it; sit=Y_structure->succ(sit))
  273.   { seg = Y_structure->key(sit);
  274.  
  275.     double cross_y = (get_slope(seg) * get_x(p) + get_yshift(seg));
  276.  
  277.     v = get_left_node(seg);
  278.     if (v==nil)
  279.     { w = New_Node(SUB,get_x(p),cross_y);
  280.       set_left_node(seg,w);
  281.      }
  282.     else
  283.     { double vx = SUB[v].xcoord();
  284.       if ( vx < get_x(p)-EPS2) 
  285.       { w = New_Node(SUB,get_x(p),cross_y);
  286.         New_Edge(SUB,v,w,seg);
  287.         set_left_node(seg,w);
  288.        }
  289.       else w = v;
  290.      }
  291.  
  292.     u = get_left_node(l);
  293.     if (u!=nil && u!=w) New_Edge(SUB,u,w,l);
  294.     set_left_node(l,w);
  295.  
  296.    }
  297.   
  298.   delete l;
  299.   delete top;
  300.   delete bot;
  301.     
  302.   Y_structure->del_item(bot_it);
  303.   Y_structure->del_item(top_it);
  304.  
  305.  }
  306.  
  307.  
  308. void Sweep_Segments(const list<segment>& L1, const list<segment>& L2, 
  309.                     GRAPH<point,int>& SUB,
  310.                     priority_queue<seq_item,sw_point>& Q,
  311.                     sortseq<sw_segment,pq_item>& S)
  312. {
  313.   sw_point    p,inter;
  314.   sw_segment  seg, l,lsit,lpred,lsucc,lpredpred;
  315.   pq_item  pqit,pxmin;
  316.   seq_item sitmin,sit,sitpred,sitsucc,sitpredpred;
  317.  
  318.   X_structure = &Q;
  319.   Y_structure = &S;
  320.  
  321.   int count=1;
  322.  
  323.   //initialization of the X-structure
  324.  
  325.   segment s;
  326.  
  327.   forall(s,L1) 
  328.    { sw_point p = new sw_POINT(s.start());
  329.      sw_point q = new sw_POINT(s.end());
  330.      seg = new sw_SEGMENT(p,q,0,count++);
  331.      Xinsert(nil,get_startpoint(seg));
  332.    }
  333.  
  334.   count = -1;
  335.  
  336.   forall(s,L2) 
  337.    { sw_point p = new sw_POINT(s.start());
  338.      sw_point q = new sw_POINT(s.end());
  339.      seg = new sw_SEGMENT(p,q,1,count--);
  340.      Xinsert(nil,get_startpoint(seg));
  341.    }
  342.  
  343.  
  344.   count = 0;
  345.  
  346.   x_sweep = -MAXINT;
  347.   y_sweep = -MAXINT;
  348.  
  349.  
  350.   while(!X_structure->empty())
  351.   {
  352.     pxmin = X_structure->find_min();
  353.     p = X_structure->inf(pxmin);
  354.  
  355.     sitmin = X_structure->key(pxmin);
  356.  
  357.     Xdelete(pxmin);
  358.  
  359.  
  360.  
  361.     if (get_kind(p) == Leftend)
  362.  
  363.     //left endpoint
  364.     { 
  365.  
  366.       l = get_seg(p); 
  367.  
  368.       x_sweep = get_x(p);
  369.       y_sweep = get_y(p);
  370.  
  371.  
  372.       if (get_x(p) == get_x(get_endpoint(l)))
  373.         vertical_segment(SUB,l);
  374.       else
  375.       {
  376.  
  377. /*
  378.       sit = Y_structure->lookup(l);
  379.       if (sit!=nil)  
  380.            error_handler(1,"plane sweep: sorry, overlapping segments");
  381. */
  382.  
  383.  
  384.       sit = Y_structure->insert(l,nil);
  385.  
  386.       Xinsert(sit,get_endpoint(l));
  387.  
  388.       sitpred = Y_structure->pred(sit);
  389.       sitsucc = Y_structure->succ(sit);
  390.  
  391.       if (sitpred != nil) 
  392.       { if ((pqit = Y_structure->inf(sitpred)) != nil)
  393.           delete Xdelete(pqit);
  394.  
  395.         lpred = Y_structure->key(sitpred);
  396.  
  397.         Y_structure->change_inf(sitpred,nil);
  398.  
  399.         if (intersection(lpred,l,inter))
  400.             Y_structure->change_inf(sitpred,Xinsert(sitpred,inter));
  401.       }
  402.  
  403.  
  404.       if (sitsucc != nil)
  405.       { lsucc = Y_structure->key(sitsucc);
  406.         if (intersection(lsucc,l,inter))
  407.            Y_structure->change_inf(sit,Xinsert(sit,inter));
  408.       }
  409.      } /* else if vertical */
  410.  
  411.     }
  412.     else if (get_kind(p) == Rightend)
  413.          //right endpoint
  414.          { 
  415.            x_sweep = get_x(p);
  416.            y_sweep = get_y(p);
  417.  
  418.            sit = sitmin;
  419.  
  420.            sitpred = Y_structure->pred(sit);
  421.            sitsucc = Y_structure->succ(sit);
  422.  
  423.            sw_segment seg = Y_structure->key(sit);
  424.  
  425.            Y_structure->del_item(sit);
  426.  
  427.            delete seg;
  428.  
  429.            if((sitpred != nil)&&(sitsucc != nil))
  430.            {
  431.              lpred = Y_structure->key(sitpred);
  432.              lsucc = Y_structure->key(sitsucc);
  433.              if (intersection(lsucc,lpred,inter))
  434.                 Y_structure->change_inf(sitpred,Xinsert(sitpred,inter));
  435.            }
  436.          }
  437.          else
  438.          /*point of intersection*/
  439.          { 
  440.            node w = New_Node(SUB,get_x(p),get_y(p));
  441.  
  442.            count++;
  443.  
  444.            /* Let L = list of all lines intersecting in p 
  445.  
  446.               we compute sit     = L.head();
  447.               and        sitpred = L.tail();
  448.  
  449.               by scanning the Y_structure in both directions 
  450.               starting at sitmin;
  451.  
  452.            */
  453.  
  454.            /* search for sitpred upwards from sitmin: */
  455.  
  456.            Y_structure->change_inf(sitmin,nil);
  457.  
  458.            sitpred = Y_structure->succ(sitmin);
  459.  
  460.  
  461.            while ((pqit=Y_structure->inf(sitpred)) != nil)
  462.            { sw_point q = X_structure->inf(pqit);
  463.              if (compare(p,q) != 0) break; 
  464.              X_structure->del_item(pqit);
  465.              Y_structure->change_inf(sitpred,nil);
  466.              sitpred = Y_structure->succ(sitpred);
  467.             }
  468.  
  469.  
  470.            /* search for sit downwards from sitmin: */
  471.  
  472.            sit = sitmin;
  473.  
  474.            seq_item sit1;
  475.            
  476.            while ((sit1=Y_structure->pred(sit)) != nil)
  477.            { pqit = Y_structure->inf(sit1);
  478.              if (pqit == nil) break;
  479.              sw_point q = X_structure->inf(pqit);
  480.              if (compare(p,q) != 0) break; 
  481.              X_structure->del_item(pqit);
  482.              Y_structure->change_inf(sit1,nil);
  483.              sit = sit1;
  484.             }
  485.  
  486.  
  487.  
  488.            // insert edges to p for all sw_segments in sit, ..., sitpred into SUB
  489.            // and set left node to w 
  490.  
  491.            lsit = Y_structure->key(sitpred);
  492.  
  493.            node v = get_left_node(lsit);
  494.            if (v!=nil) New_Edge(SUB,v,w,lsit);
  495.            set_left_node(lsit,w);
  496.  
  497.            for(sit1=sit; sit1!=sitpred; sit1 = Y_structure->succ(sit1))
  498.            { lsit = Y_structure->key(sit1);
  499.  
  500.              v = get_left_node(lsit);
  501.              if (v!=nil) New_Edge(SUB,v,w,lsit);
  502.              set_left_node(lsit,w);
  503.             }
  504.  
  505.            lsit = Y_structure->key(sit);
  506.            lpred=Y_structure->key(sitpred);
  507.            sitpredpred = Y_structure->pred(sit);
  508.            sitsucc=Y_structure->succ(sitpred);
  509.  
  510.  
  511.            if (sitpredpred != nil)
  512.             { 
  513.               lpredpred=Y_structure->key(sitpredpred);
  514.  
  515.               if ((pqit = Y_structure->inf(sitpredpred)) != nil)
  516.                 delete Xdelete(pqit);
  517.  
  518.               Y_structure->change_inf(sitpredpred,nil);
  519.  
  520.  
  521.               if (intersection(lpred,lpredpred,inter))
  522.                 Y_structure->change_inf(sitpredpred,
  523.                                        Xinsert(sitpredpred,inter));
  524.              }
  525.  
  526.  
  527.            if (sitsucc != nil)
  528.             {
  529.               lsucc=Y_structure->key(sitsucc);
  530.  
  531.               if ((pqit = Y_structure->inf(sitpred)) != nil)
  532.                 delete Xdelete(pqit);
  533.                  
  534.               Y_structure->change_inf(sitpred,nil);
  535.  
  536.               if (intersection(lsucc,lsit,inter))
  537.                   Y_structure->change_inf(sit,Xinsert(sit,inter));
  538.              }
  539.  
  540.  
  541. // reverse the subsequence sit, ... ,sitpred  in the Y_structure
  542.  
  543.            x_sweep = get_x(p);
  544.            y_sweep = get_y(p);
  545.  
  546.            Y_structure->reverse_items(sit,sitpred);
  547.  
  548.           delete p;
  549.  
  550.          } // intersection
  551.  
  552.   }
  553.  
  554.   X_structure->clear();
  555.  
  556. }
  557.  
  558.  
  559.  
  560. #if !defined(__TEMPLATE_ARGS_AS_BASE__)
  561. declare3(_sortseq,sw_segment,pq_item,skiplist)
  562. declare3(_sortseq,sw_segment,pq_item,ab_tree)
  563. #endif
  564.  
  565. main()
  566. { 
  567.   int N = read_int("N = ");
  568.  
  569.   list<segment> seglist1,seglist2;
  570.   
  571.   while (N--) 
  572.   { double x1 = random(-1000,-100);
  573.     double y1 = random(-1000,1000);
  574.     double x2 = random(100,1000);
  575.     double y2 = random(-1000,1000);
  576.     seglist1.append(segment(x1,y1,x2,y2));
  577.    }
  578.  
  579.   GRAPH<point,int> SUB;
  580.  
  581.   sortseq<sw_segment,pq_item>           rst_seq;
  582.  
  583. #if defined(__TEMPLATE_ARGS_AS_BASE__)
  584.   _sortseq<sw_segment,pq_item,ab_tree>   abt_seq;
  585.   _sortseq<sw_segment,pq_item,skiplist>  skip_seq;
  586. #else
  587.   _sortseq(sw_segment,pq_item,ab_tree)   abt_seq;
  588.   _sortseq(sw_segment,pq_item,skiplist)  skip_seq;
  589. #endif
  590.  
  591.   priority_queue<seq_item,sw_point>     Q;
  592.  
  593.  
  594.   float T;
  595.  
  596.   T = used_time();
  597.   cout << "SWEEP_SEGMENTS: ";
  598.   cout.flush(); 
  599.   SWEEP_SEGMENTS(seglist1,seglist2,SUB);
  600.   cout<< string(" # = %d time = %6.2f sec",SUB.number_of_nodes(), used_time(T));
  601.   newline;
  602.  
  603.   SUB.clear();
  604.  
  605.   T = used_time();
  606.   cout << "rs_tree:        ";
  607.   cout.flush(); 
  608.   Sweep_Segments(seglist1,seglist2,SUB,Q,rst_seq);
  609.   cout<< string(" # = %d time = %6.2f sec",SUB.number_of_nodes(), used_time(T));
  610.   newline;
  611.  
  612.   SUB.clear();
  613.  
  614.   T = used_time();
  615.   cout << "ab_tree:        ";
  616.   cout.flush(); 
  617.   Sweep_Segments(seglist1,seglist2,SUB,Q,abt_seq);
  618.   cout<< string(" # = %d time = %6.2f sec",SUB.number_of_nodes(), used_time(T));
  619.   newline;
  620.  
  621.   SUB.clear();
  622.  
  623.   T = used_time();
  624.   cout << "skiplist:       ";
  625.   cout.flush(); 
  626.   Sweep_Segments(seglist1,seglist2,SUB,Q,skip_seq);
  627.   cout<< string(" # = %d time = %6.2f sec",SUB.number_of_nodes(), used_time(T));
  628.   newline;
  629.  
  630. return 0;
  631.  
  632. }
  633.